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'"^ . Abstract. We report the results of three dimensional molecular dynamic simulations 

C I on the angle of repose of a sandpile formed by pouring mono-sized cohesionless spherical 

grains into a granular Hele-Shaw cell. In particular, we are interested to investigate 
the effects of those variables which may impact significantly on pattern formation of 
granular mixtures in Hele-Shaw cells. The results indicate that the frictional forces 
^ ' influence remarkably the formation of pile on the grain level. Furthermore, We see that 

CN I increasing grain insertion rate decreases slightly the angle of repose. We also find that 

^T" ' in accordance with experimental results, the cell thickness is another significant factor 

p^ . and the angle of repose decays exponentially by increasing the cell thickness. It is 

shown that this effect can be interpreted as a cross-over from two to three dimensions. 
In fact, using grains with different sizes shows that the behavior of the angle of repose 
C^^ , when both size and cell thickness are varied is controlled by a scaled function of the 

^— ^ ' ratio of these two variables. 
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1. introduction 

The formation of a sandpile is an everyday life experience which may occur in several 
ways, including pouring |lj, discharging [21 [3], failing [HE], tilting and rotating [6l[7]. 
This gravity- driven physical event is relevant to many practical applications and 
laboratory processes such as avalanches [H [9l [lO] and pattern formation in granular 
mixtures [HI [12]. It is also related to the paradigm of self organized criticality [13] which 
has itself a wide realm of applicability to a diverse range of physical phenomena [Ti] . 

The angle of repose of a sandpile is a simple concept that characterizes the behaviour 
of granular materials on macroscopic scale. A sandpile has an inclined free surface that 
does not flow. When grains are added to the sandpile the slope of this free surface 
increases until, it exceeds a threshold value. At this point, the pile does not support 
new grains and releases many grains in an avalanche reducing the slope to the angle 
of repose, 6r- One of the most important questions in the study of piling of grains 
is how to determine the relationship between the angle of repose and the relevant 
variables, i.e. material properties, boundary conditions, and the model interactions 
at the grain level. Experiments have shown that 6r depends on the shape, density, and 
size distribution of the grains, as well as humidity, formation history and geometrical 
boundary conditions [Il[2l[3]. 

Granular materials are ensembles of discrete particles, and as for any macroscopic 
system, the total number of modes are huge. What we see on macroscopic scale, hence, 
is a manifestation of the collective behaviour of a large assembly of macroscopic grains 
that interact with each other (and maybe with the container walls), through collisions 
and friction. In this regard, the molecular dynamics (MD) simulation technique provides 
perhaps the most realistic approach for modeling the properties of granular materials. 
The present work is focused on such a numerical study of the angle of repose of a 
heap formed by pouring dry, mono-sized spherical grains into a granular Hele-shaw cell. 
Some of the experimental aspects of the phenomenon have been reported by Grasselli 
and Herrmann [l]. However, as far as we know, there has been no three-dimensional 
MD study of the angle of repose in Hele-Shaw cell, except the one by Zhuo et al [15], 
which studied the formation of sandpiles resulting from discharging materials by using 
a different model interaction. Regarding to the classification made by Grasselli and 
Herrmann [1] , what they have estimated is the outflow angle of the heap which is less 
than the heap angle, the subject of current research. 

In this paper, we report a detailed numerical study of the formation of a sandpile 
and the dependence of its angle of repose on the relevant variables. We are specifically 
interested to determine how 6r changes with variables like the size of particles, the 
rate of mass injection, the density of grains, the thickness of the container, and wall- 
particle and particle-particle static friction coefficients. These are the key parameters 
which may impact significantly on pattern formation of granular mixtures in Hele-Shaw 
cells [ini [17]. To answer such questions, we performed extensive MD simulations in 
three dimension on model systems by using a new MD scheme developed by Silbert et 



The angle of repose of spherical grains... 



Particle 




Particle j 



Figure 1. Two-dimensional illustration of two grains i and j in contact and position 
vectors r^jFj, respectively, with compression and overlap Sij. 



al [18], which is described briefly in the next section. The results of our simulations are 
presented in section 3, followed by the main conclusions at section 4. 

2. SIMULATION METHODOLOGY 

As mentioned earlier, our simulations are based on a Distinct Element Method (DEM) 
scheme, originally proposed by Silbert et al. The method has been discussed elsewhere 
in details p^. For the sake of clarity, we review the most important aspects of it here. 
According to this model, the spheres interact on contact through a linear spring dashpot 
or Hertzian interactions in the normal and tangential directions to their lines of centers. 
In a gravitational fleld g, the translational and rotational motions of grain z in a system 
at time t, caused by its interactions with neighboring grains or walls, can be described 
by Newton second law, in term of the total force and torque as the following equations: 

Fr="!g+E (%■+%-) (1) 

Where mi,rij = Fj — r^ are respectively, the mass of grain i and the relative distance 
between grains i and j. For two contacting grains i, j at positions r^, Tj with velocities 
Vj,Vj and angular velocities uii,ujj the force on grain i is computed as follows: the 
normal compression 6ij is (FigHj) 

kj = d- r^^j (3) 

and relative normal velocity v„^ . and relative tangential velocity v^. , are given by 



Vr 



'i ,j \ ^i3 



•ni,i)n.j (4) 



vt,,, = Vij - v„^^^ - -(^i + ^i)rij (5) 
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Where rijj- = ^^ with rij = |rj j| and Vjj- = Vj — Vj. The rate of change of elastic 
tangential displacement u^. ,, set to zero at the initial of a contact, is given by [18] 

dUf. . (u+. ..Vj ~)Ti j , , 

The second term in Eq.(6) arises from the rigid body rotation around the contact 
point and insures that u^- . always lies in the local tangent plan of contact. In Eqs. (1) 
and (2), normal and tangential forces acting on grain i are given by 

Fn,,, = f{^){kJi,jnij - -fnmeffVn^J (7) 

and 

Ft,,, = fi^)iktut^^ - -itm.ffWtJ (8) 

Where k^^t and 7„.t are elastic and viscoelastic constants respectively and ?7ie// = 
mimj/{mi + rrij) is the effective mass of spheres with masses nii and ruj. The 
corresponding contact force on grain j is simply given by Newtons third law, i.e., 
Fj J- = — Fj j . For spheres of equal mass m, as is the case in our system, me// = m/2; 
f{x = 1) for the linear spring dashpot (Hookian) model, and f{x) = ^Jx for Hertzian 
contacts with viscoelastic damping between spheres [T8] . 

Static friction is implemented by keeping track of the elastic shear displacement 
throughout the lifetime of a contact. The static yield criterion, characterized by a 
local particle friction coefficient /i, is modeled by truncating the magnitude of u^^ as 
necessary to satisfy |Fi. J < |/iF„. J. Thus the contact surfaces are treated as sticking 
when |Fi- .| < |/iF„. ,|, and as slipping when the yield criterion is satisfied p^ [20]. 



3. RESULTS AND DISCUSSION 

In this section we present the results of our extensive molecular dynamics (MD) 
simulations in three dimensions on model systems of N mono- disperse, cohesionless and 
inelastic spheres of diameter d and mass m. The system is constrained by a rectangular 
box with fixed rough walls and boundary conditions and free top surface, as in FiglJ] 
A simulation was started with the random generation of spheres without overlaps from 
top and left corner of container, followed by a gravitational settling process to form 
a stable heap. The results are given in non-dimensional quantities by defining the 
following normalization parameters: distance, time, velocity, forces, elastic constants, 
and stress are, respectively measured in units of d, to = yd/ g, vq = \/dg, Fq = mg, and 
fco = mg/d. All data was taken after the system had reached the steady state. Because 
of the complexity of the model, there are a wide range of parameters that affect the 
results of computation. However, we usually investigate the effect of a single variable 
varying in a certain range while other variables are fixed to their base values as listed 
in table [H 

All the cases were simulated in three dimensions using a molecular dynamics code 
for granular materials LAMMPS [HI [21]. The equations of motion for the translational 
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Figure 2. Geometry of cell used in simulations. 



and rotational degrees of freedom are integrated with either a thirdorder Gear predictor- 
corrector or velocity- Verlet scheme [22] . The angle of repose could then be determined 
from the surface profile of the heap. 



3.1. The evolution of the pile 

We start from an empty vertical cell composed of two parallel plates separated by an 
spacer with a variable width w. A granular heap is formed in the cell by pouring mono- 
sized spherical grains on the bottom plate. These grains are released from a small box 
located on the top-left of the cell, as shown in FiglH The rate of material insertion R 
can be varied by changing the box size. A pile with well defined shape starts to form as 
the number of the added grains growth. The angle of repose could then be determined 
from the surface profile of the heap. 

The formation of a granular heap at the macroscopic level is a complicated 
phenomenon arises as the result of surface avalanches and also a kink mechanism |23j . 
At the grain level, it is a consequence of energy dissipation. In addition to frictional 
forces, the particle-particle inelastic collisions consumes the translational and rotational 
energy of a falling grain. Eventually, the grain will stick to a position where its energy 
is less than the minimum value required to overcome the potential barrier created by 
frictional forces ^24j. At the earlier stage of the pile formation, the number of grains are 
small. Therefore, the falling grains can move further and the angle of repose is less than 
the one of a larger pile. By adding more grains to the pile the angle of repose growth to 
its asymptotic value, where the energy gain in the gravitational field is balanced by the 
dissipating effects. In Figl3]the evolution of the sandpile or equivalently, the evolution of 
6r has been plotted for two values of fi^ against A^, the total number of released grains. 
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Table 1. Basic computational parameters 



Parameters system conditions 

Number of Particles (N) 40,000 

Particle Size (d) Ido 

Particle density (p) 2.4(1710/(1^) 

Particle friction coefficient (/ip) 0.5 

Wall friction coefficient (fiw) 0.5 

Particle normal stiffness coefficient (fc„) 2 x lO'^(fco) 
Particle tangential stiffness coefficient (kt) jkn 

Particle normal damping coefficient (7„) 50/(io) 

Particle tangential damping coefficient (7^) ^{^0 ) 

Wall normal stiffness coefficient (fc„) 2 x lO'^(fco) 

wall tangential stiffness coefficient (fcj) |/c„ 

wall normal damping coefficient {"/„) 50/(io) 

Particle normal damping coefficient (74) 50/(io) 

Time step increment 2 x 10^'^ 



As seen from the figure, the angle of repose grows hnearly for smaller values of A^ and 
then evolves quickly into a steady state, characterized by an asymptotic value of 6r at 
larger values of A^. For the set of parameters used in FiglS], the angle of repose is almost 
constant for A^ > 25, 000. For the rest of the work we have used A^ = 40, 000 grains in 
all simulations to make sure that the angle of repose has reached its asymptotic values 
as shown in Figj3l 

We also note that at earlier stage of heap formation the difference between 6{N)^s 
is not significant for different values of fj,.^, as seen in FigjSl This maybe explains 
the existence of the so called dead-zone [12] in stratification of bi-dispersed mixtures. 
It is known that the difference between the angle of repose of two species is the key 
factor which together with the difference between the grain sizes determines the type of 
patterns formed in the cell [TT]. If the difference between 6r is negligible, e.g., when the 
number of grains is small, the only factor that affects pattern formation is the size. As 
such, for small values of A^, the grain segregate according to their size and a dead- zone 
forms. 

3.2. Effect of mass flow rate 

A simple way of increasing the rate of mass flow into the cell is increasing the density of 
each grain, while its diameter is fixed. From FigJUone can see that increasing the mass 
of individual grains causes a slight fluctuation, at most 0.3° around the mean values 
in both linear and Hertzian model. In fact, there is no hint of a dependence of 6r on 
the value of grain's mass in our simulations in the range of parameters we used. This 
conclusion is in agreement with some experimental observations [TT] . 

A more interesting way of increasing mass flow rate is to increase R, the number of 
grains falling down on the heap at each time step. It has been found that the wavelength 
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Figure 3. Angle of repose as a function of number of grains from Hookin model. 



of the layers in granular stratification is an increasing function of this quantity 
Moreover, above a critical flow rate Re, the stratification pattern abruptly disappear [T7|. 
Our simulations shown in Figj5] indicate that, increasing the rate of grain insertion R, 
decreases the angle of repose slightly. The change in 6r is small, about a degree after 
increasing R three times. Such a slight change in the value of angle of repose when the 
rate of grain insertion was varied has been reported by Grasselli and Herrmann [T] . 

Comparing the two ways of increasing R, we conclude that in the formation of a 
sandplie, what does really matter is not the mass itself, but the number of grains that 
hit the pile in a given time. This effect may be attributed to the occurrence of larger 
avalanches when the rate of grains falling upon top of the heap becomes larger. 



3.3. Effect of cell thickness 

The effect of variation in the size of the gap between walls container on the piling of 
grains has been the subject of both experimental and numerical studies with the general 
conclusion that increasing the cell thickness induces an exponential decay in the angle 
of repose [1] 



6r(w) = 6oc(l - ae 



-w/l\ 



(9) 



where a is a constant depending on the grain properties and I is a characteristic length 
representing the scale over which the walls affect the piling of grains. The parameter w 
is also a determining factor which affects drastically pattern formation of bi-dispersed 
granular mixtures in Hele-Shaw cells. A nonlinear decrease in the wavelength of the 
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Figure 4. The angle of repose against density from Hookian and Hertzian Models 



on 



o 
a. 

a 

O 



28 



27 - 



26 



25 - 



24 - 



(.01 



. 














. 








D 


Hookian Model, w=8 


- 








A 


Hertzian Model, w=8 


- 








O 


Hookian Model, w=16 


- 








o 


Hertzian Model, w=16 


- 


D 












- 


A 


D 


D 








. 




A 


A 


R 






. 










S 




- 












^ 


- 


<§> 


9 

1 , , 


8 

1 , 


8 

, , 1 , 


O 


o 
o 



0.015 0.02 0.025 0.03 

Particle insertion/time step 



0.035 



Figure 5. Angle of repose as a function of grain insertion from Hookian and Hertzian 
models for two different container thickness. 
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Figure 6. The angle of repose against cell thickness simulated with Hookian and 
Hertzian Models. 



stratification pattern with cell thickness w has been reported [16]. And similar to Re, 
there is a critical cell thickness Wc above which the layered pattern abruptly vanishes |17j . 

To quantify the behaviour of the angle of repose when the cell thickness is varied, 
simulations were performed using different cell thickness w. The results presented in 
Figl6] shows that by increasing the cell thickness, the angle of repose, 9r{w), decreases 
and it eventually settles down to an asymptotic value ^oo when the cell thickness becomes 
very large compared to the size of grains. A logarithmic plot of (^oo — 0r{w))/6oo 
versus w shows that in accordance with experimental observations, the trend follows an 
exponential decay (FigJT]). We also observe that for both Hookian and Hertzian models, 
have the same slope, which means that the characteristic lengths does not depend on 
the model interaction. Since MD simulations on angle of repose of granular heap formed 
by discharging materials have produced the same results [15], one may conclude that 
the observed phenomenon is universal and does not depend in the formation history of 
the pile. 

This phenomenon can be explained based on a dimensionality analysis. The walls of 
the cell not only consume the kinetic energy of grains via friction and inelastic collisions, 
but they confine them to move in a limited space. Indeed, by changing the cell thickness, 
while keeping the size of grains fixed, one can even change the dimensionality of the 
space accessible to the grains. Investigation of the average number of contacts per 
grain, z, in this round of simulations demonstrates the validity of this image. In FigJS] 
the variation of 1 versus w has been depicted for the interior grains, i.e. grains which 
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Figure 8. The average number of contacts for interior grains versus cell thickness 
from Hookian model. 
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are not supported by the walls. The fraction of interior grains decreases with w/d, but 
even at a small value of cell thickness as w/d = 1.5, there are still a non-negligible 
fraction of such grains (about %8.8 of total grains) in the cell. We also observe that 
for w > 2d, the average number of contacts becomes less than 4, which is the minimal 
average coordination number required to obtain a static packing of frictional spheres in 
three dimension [251 126] indicating a cross-over from a three dimensional regime to a 
two dimensional one. At w = 1.5d, the mean number of compact is 3.3, very close to the 
minimal value z = 3 needed for a stable packing of frictional spheres in two dimension. 
On the other hand, for large values of w/d, where the angle of repose reaches to its bulk 
value, the value of z quickly reaches to the asymptotic value z = 5.2. Of course, this 
limiting value depends on the microscopic parameters of the granular materials. 

3.4- Effect of particle size 

Another important and interesting parameter which we expect to affect the piling of 
grains is the grain diameter d itself. In general, one expects that the angle of repose of a 
three dimensional sandpile does not depend on the grain size, since a pile of smaller 
grains transforms into a pile of larger grains via an isotropic re-scaling of the pile 
coordinations [16]. Based on this idea, we expect that in granular Hele-Shaw cells, 
the angle of repose of a pile to be a function of w/d, because the confinement of the pile 
between two vertical plates, remove the spatial symmetry in the perpendicular direction. 
As such, in an isotropic re-scaling of a sandpile in a Hele-Shaw cell, the gap between 
walls should be re-scaled with grain size too. 

In order to examine this idea, we performed a series of simulations using the base 
parameter values as listed in table. I, but with different different grain size d and cell 
thickness w, using linear model. The results are sketched in Figj9l where the behaviour 
of 9r has been plotted against w/d, for several values of grain diameters. As anticipated, 
all the curves almost collapse into one curve showing that in general, the angle of repose 
of a sandpile in a Hele-Shaw cell does depend on the size of grains. In the range of 
parameter values we used here, the value of 6r decays quickly to its asymptotic value 
^oo for w/d > 20. This is in agreement with our anticipation that in three dimension 
the angle of repose does not depend on the size of grains. The data collapse in Figj9] 
implies that the parameters a and 6oo in Eq.(ll) do not depend on the grain size. 
Furthermore, the value of characteristic length should be proportional to the grain size 
d. In the experiments conducted by Grasselli and Herrmann, the value of 6oo and a turn 
out to be size-independent too. Meanwhile, they did not find any evidence for a linear 
dependence of the characteristic length on grain diameter. In fact, the characteristic 
length I in their observation was the same for all spherical grain sizes. This unexpected 
result could be due to the slightly variation in the size of grains (polydispersity of the 
granular materials) they used [1]. It would also be a result of the cohesive force between 
grains resulted from capillary force for humidified grains or from van der Waals force 
for dry grains less than 100 fim [15]. 
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Figure 9. Angle of repose as a function of cell thickness scaled with particle size in 
Hookian and Hertzian models. 



On the other hand, MD simulations on heaps formed by discharging grains indicate 
that 6r depends on particle size such that, for a given w/d, larger grain have smaller 
angle of repose |15] . This is in contrary to our results described above. But as mentioned 
earlier, not only the model interactions, but the formation histories are different in these 
two simulations. The importance of the latter factor is a well-known experimental fact. 
Moreover, there are evidences that in the absence of cohesive force between grains or 
size-dependent sliding coefficient [27] , the size-dependent angle of repose could arise as 
a result of introducing the rolling friction as Zhuo et al did so [15] . 



3.5. Effects of frictional forces 

The friction between grains exhausts the translational and rotational energy of particles. 
It also affects significantly the mechanical stability of each contact in the sandpile. 
One expects, therefore, that the particle-particle friction coefficient /ip be one of the 
key parameters in determining the angle of repose of the pile. To investigate how 
this variable change the behaviour of Or, we performed computer simulations for the 
spherical particles of fixed size with different /ip using both Hertzian and Hookian 
models. The other values are listed in table I. The results have been presented on 
FigHni showing that for both models, the angle of repose increases non-linearly with 
fip. This effect is in accordance with MD simulations of 6'^ when the pile forms from 
discharging materials [3l |15] . 

We have also examined the effect of variation of the wall-particle friction coefficient 
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Figure 10. Angle of repose as a function of friction coefficient between grains from 
Hookian and Hertzian models. 



/i^ on Or. This is another important parameter which should affect the angle of repose, 
as wall friction gives a torque resistance to the rotational motion of grains. The results 
of simulation of the wall friction coefficient are shown in FigJTTl As indicated in the 
figure, increasing the wall friction can significantly increase the angle of repose. Similar 
trends have also been observed for the grains of different sizes. 

In Fig {121 we have compared the behaviour of 9r versus fiw, for two different kn 
values. It is seen that increasing material stiffness the angle of repose gets smaller, 
but the change is not very significant. The decrease in 6r as the result of increasing 
kn was expected. But the important fact is that the change in Or after increasing the 
material stiffness coefficients by an order of magnitude is only 2° — 3°. This indicates 
that parameters like kn does not play a very crucial rule in formation of sandpiles. 



4. CONCLUSIONS 

We study the dependence of the angle of repose of a sandpile resulting from pouring 
granular materials in Hele-Shaw cell on boundary conditions, particle characteristics, 
and the other parameters used in the model interaction. We saw that the most 
important boundary condition is the cell thickness. Increasing this parameter induces 
an exponential decay in angle of repose, in agreement with experimental findings. The 
same trend has been observed in MD simulation of discharging materials, showing that 
the effect is rather universal and independent of the formation history of the pile. For 
cells with a large gap, the angle of repose relaxes to an asymptotic value which should 
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Figure 12. Comparison of the angle of repose as a function of friction coefficient 
between particles and walls from Hookian model for two different elastic coefficients. 
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be correspond to the angle of repose of an isotropic three dimensional heap. 

At the grain level, the frictional forces are among the most important factors. We 
observed that increasing particle-particle friction and wall-particle friction coefficient 
increase the angle of repose non-linearly. But the most interesting microscopic 
characteristic of the single grains maybe is the grain diameter d. We found that for 
a given cell thickness, 9r increases with grain size, but if the wall thickness scales with 
grain diameter, the angle of repose remain unchanged. This also implies that the angle 
of repose in three dimensional limit does not depend on the size of grains. It should be 
emphasized that these observations are just valid for the case of cohesionless, mono-sized 
grains. As such, they might be in contrary with the result of real experiments, where 
these conditions are not satisfied or the formation history is different. 
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